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Electron dispersion forces play a crucial role in determining the structure and properties of 
biomolecules, molecular crystals and many other systems. However, an accurate description of 
dispersion is highly challenging, with the most widely used electronic structure technique, density 
functional theory (DFT), failing to describe them with standard approximations. Therefore, appli- 
cations of DFT to systems where dispersion is important have traditionally been of questionable 
accuracy. However, the last decade has seen a surge of enthusiasm in the DFT community to tackle 
this problem and in so-doing to extend the applicability of DFT-based methods. Here we discuss, 
classify, and evaluate some of the promising schemes to emerge in recent years. A brief perspective 
on the outstanding issues that remain to be resolved and some directions for future research are also 
provided. 
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I. INTRODUCTION 

The theoretical description of matter as well as of 
many chemical, physical, and biological processes re- 
quires accurate methods for the description of atomic 
and molecular-scale interactions. Whilst there are many 
quantum mechanical approaches, in the past few decades 
Kohn-Sham density functional theory (DFT) 1 ' 2 has es- 
tablished itself as the theoretical method of choice for this 
task, undergoing a meteoric rise in large parts of physics, 
chemistry, and materials science. The rise of DFT and 
its uptake in academia and industry has been widely 
discussed, and was perhaps illustrated most clearly in 
Burke's recent Spotlight article on DFT. 3 

Although DFT is in principle exact, in practice ap- 
proximations must be made for how electrons interact 
with each other. These interactions are approximated 
with so-called exchange-correlation (XC) functionals and 
much of the success of DFT stems from the fact that 
XC functionals with very simple forms often yield accu- 
rate results. However there are situations where the ap- 
proximate form of the XC functional leads to problems. 
One prominent example is the inability of "standard" 
XC functionals to describe long-range electron correla- 
tions, otherwise known as electron dispersion forces; by 
standard XC functionals we mean the local density ap- 
proximation (LDA), generalized gradient approximation 
(GGA) functionals or the hybrid XC functionals. The 
"lack" of dispersion forces - often colloquially referred to 
as van der Waals (vdW) forces - is one of the most signif- 
icant problems with modern DFT and, as such, the quest 
for DFT-based methods which accurately account for dis- 



persion is becoming one of the hottest topics in compu- 
tational chemistry, physics, and materials science. Fig. 1 
underlines this point, where it can be seen that over 800 
dispersion-based DFT studies were reported 2011 com- 
pared to fewer than 80 in the whole of the 1990s. 
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FIG. 1. The number of dispersion corrected DFT studies has 
increased considerably in recent years. The total number of 
studies performed is difficult to establish precisely, however, 
an estimate is made here by illustrating the number of papers 
that cite at least one of 16 seminal works in the field (Refs. 
4-19). (Data from Web of Knowledge, July 2012 for the years 
1991-2011) 

Dispersion can be viewed as an attractive interaction 
originating from the response of electrons in one region 
to instantaneous charge density fluctuations in another. 
The leading term of such an interaction is instantaneous 
dipole-induced dipole which gives rise to the well known 
— 1/r 6 decay of the interaction energy with interatomic 
separation r. Standard XC functionals don't describe dis- 
persion because: (a) instantaneous density fluctuations 



are not considered; and (b) they are "short-sighted" in 
that they consider only local properties to calculate the 
XC energy. The consequence for two noble gas atoms, for 
example, is that these functionals give binding or repul- 
sion only when there is an overlap of the electron densities 
of the two atoms. Since the overlap decays exponentially 
with the interatomic separation, so too does any binding. 
We show this in Fig. 2 for one of the most widely used 
GGAs - the so-called Perdew-Burke-Ernzerhof (PBE) 
functional 20 - for a binding curve between two Kr atoms. 
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FIG. 2. Binding curves for the Kr dimer obtained with the 
PBE exchange-correlation functional and an accurate model 
potential. 21 Dispersion originates from fluctuations in the 
electron density which polarize different atoms and molecules, 
as shown schematically in the lower right diagram. The in- 
teraction then exhibits the well-known — 1/r 6 decay. This 
— 1/r 6 decay is not reproduced with PBE (or other semi-local 
functionals) which instead gives an exponential decay for the 
interaction because PBE relies on the overlap of density to 
obtain the interaction. 

The binding of noble gases is a textbook dispersion 
bonded system. However, it has become increasingly 
apparent that dispersion can contribute significantly to 
the binding of many other types of materials, such as 
biomolecules, adsorbates, liquids, and solids. Fig. 3 illus- 
trates a more "real world" example, where an accurate 
description of dispersion is critical. The figure reports 
binding energies obtained from PBE and an accurate ref- 
erence method for DNA base pairs in two different config- 
urations. In the first configuration the binding between 
the base pairs is dominated by hydrogen bonding. Hy- 
drogen bonding is governed mainly by electrostatics and 
a standard functional such as PBE predicts reasonable 
hydrogen bond strengths and as a result the stability of 
the dimer is quite close (within 15%) to the reference 
value. 22 However, in the other "stacked" configuration 
the binding is dominated by dispersion forces and for 
this isomer PBE is hopeless, hardly predicting any bind- 
ing between the base pairs at all. This huge variability 
in performance is far from ideal and since the stacked ar- 
rangement of base pairs is a common structural motif in 
DNA, the result suggests that DNA simulated with PBE 
would not be stable. 

As we will see throughout this article there are many 
other examples of the importance of dispersion to the 
bonding of materials and in recent years a plethora of 
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FIG. 3. Two binding configurations of the DNA base pairs 
adenine and thymine. A hydrogen bonded structure is shown 
on the left (hydrogen bonds indicated by red dots) and a 
"stacked" geometry on the right. For the hydrogen bonded 
configuration a standard XC functional such as PBE gives a 
binding energy that is in rather good agreement with the ref- 
erence value (lower part of the graph, data in kJ/mol). 22 In 
contrast, PBE fails to describe the binding of the "stacked" 
configuration where dispersion interactions contribute signif- 
icantly to the intermolecular binding. 



schemes has been proposed to treat dispersion within 
DFT. Here we will discuss some of the main approaches 
developed and in the process attempt to provide a use- 
ful classification of them. We also highlight some of the 
obstacles that must be overcome before improved DFT- 
based methods for including dispersion are available. By 
its nature, this article is a limited personal view that can- 
not cover every development in this thriving field. For 
the most part we have tried to keep the overview simple, 
and have aimed it primarily at newcomers to the field, 
although we do go into more depth at the later stages. 
For more detailed discussions of some of the methods 
shown here the interested reader should consult the re- 
views of Grimme, 23 Tkatchenko et al. 2A Johnson et al. 25 
or others. 26 29 Some relevant developments are also dis- 
cussed in the reviews of Sherrill 30 and Riley et al. 31 which 
are more focused on post Hartree-Fock (HF) methods. 



II. A CLASSIFICATION OF THE COMMON 
DFT-BASED DISPERSION METHODS 

Many DFT-based dispersion techniques have been de- 
veloped and rather than simply listing them all, it is use- 
ful to try to classify them. A natural way to do this is to 
consider the level of approximation each method makes in 
obtaining the long range dispersion interactions, that is, 
the interactions between well separated fragments where 
dispersion is clearly defined. In doing this it turns out 
that groups of methods which exploit similar approxi- 
mations emerge. Here, we simply identify these group- 
ings and rank them from the most approximate to the 
more sophisticated approaches in a manner that loosely 
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resembles the well-known "Jacob's ladder" of functionals 
introduced by Perdew. 32 Therefore, by analogy, we in- 
troduce a "stairway to heaven" for long range dispersion 
interactions and place each group of dispersion correction 
schemes on a different step of the stairway. In complete 
analogy to the ladder, when climbing the stairway pro- 
gressively higher overall accuracy can be expected until 
exact results, and thus heaven, are reached. We stress 
the point "overall accuracy" because, as with the lad- 
der, climbing the stairway does not necessarily mean 
higher accuracy for every particular problem but rather 
a smaller probability to fail, i.e., a statistical improve- 
ment in performance. 33 A schematic illustration of the 
stairway is shown in Fig. 4 and in the following sections 
each step on the stairway is discussed. 




persion plays a major role such as graphite or noble gases 
on metals. However, the results with LDA for dispersion 
bonded systems have limited and inconsistent accuracy 
and the asymptotic form of the interaction is incorrect. 
More promising approaches on the ground level of our 
stairway are density functionals specifically fitted to re- 
produce weak interactions around minima as well as spe- 
cially adapted pseudopotentials. 

The "Minnesota functionals" 34 are an example of a 
new breed of functionals that are fitted to a dataset that 
includes binding energies of dispersion bonded dimers, 
amongst other properties. Although these functionals 
can describe binding accurately at separations around 
minima, they suffer from the same incorrect asymptotics 
as the LDA does. On the plus side, the reference data 
used in their construction also contains properties other 
than weak interactions (e.g. reaction barriers), so that 
they can be rather accurate for general chemistry prob- 
lems. This is a clear advantage over "proper" dispersion 
correction schemes which often utilise GGA functionals 
that can be less accurate for such problems. 

For electronic structure codes that make use of pseu- 
dopotentials, dispersion can also be modelled by adding a 
specially constructed pseudopotential projector. Within 
this class are the dispersion corrected atom-centered 
potentials (DCACP) 11 and the local atomic potentials 
(LAP) methods. 35 These approaches have shown a lot 
of promise, 36 however, effort is required to carefully fit 
the potentials for each element and they are not easily 
transferable to all-electron methods. 



FIG. 4. In analogy with the Jacob's ladder classification of 
functionals the "stairway to heaven" is used here to clas- 
sify and group DFT-based dispersion correction schemes. At 
ground level are methods which don't describe the long range 
asymptotics. Simple Cq correction schemes sit on the first 
step, on the second step are approaches that utilise environ- 
ment dependent Cq corrections. The long range density func- 
tionals sit on step 3 and on step 4 and above are approaches 
which go beyond pairwise additive determinations of disper- 
sion. Upon climbing each step of the stairway the level of ap- 
proximation is reduced and the overall accuracy is expected 
to increase. 



A. Ground — Binding with incorrect asymptotics 

First, the ground is occupied by methods that sim- 
ply don't describe the long range asymptotics. These 
approaches give incorrect shapes of binding curves and 
underestimate the binding of well separated molecules. 
Some of these methods are, nevertheless, being used for 
weakly bonded systems. Although this might seem odd, 
this somewhat misguided approach is used because some 
standard DFT functionals bind dispersion bonded sys- 
tems at short separations. A prominent example is the 
LDA which has been used to study systems where dis- 



B. Step one — Simple Cq corrections 

The basic requirement for any DFT-based dispersion 
scheme should be that it yields reasonable — 1/r 6 asymp- 
totic behavior for the interaction of particles in the gas 
phase, where r is the distance between the particles. A 
simple approach for achieving this is to add an additional 
energy term which accounts for the missing long range 
attraction. The total energy then reads 

^tot = ^DFT + ^disp j (1) 

where £?dft is the DFT total energy computed with a 
given XC functional. The dispersion interaction is given 

by 

£disp = -£ C 6 AB Ai B , (2) 
A,B 

where the dispersion coefficients C$ B depend on the el- 
emental pairs A and B. Within this approach dispersion 
is assumed to be pairwise additive and can therefore be 
calculated as a sum over all pairs of atoms A and B. 
Methods on step 1 of the stairway use coefficients that 
are tabulated, isotropic (i.e. direction independent) and 
constant, and these methods are generally termed "DFT- 
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D". 

Mostly because of the simplicity and low computa- 
tional cost this pairwise Cq/t 6 correction scheme is 
widely used. Nonetheless, it has at least four clear short- 
comings which limit the accuracy one can achieve with it. 
First, the Cq/t 6 dependence represents only the leading 
term of the correction and neglects both many-body dis- 
persion effects 37 and faster decaying terms such as the 
Cs/r 8 or Cio/r 10 interactions. Second, it is not clear 
where one should obtain the Cq coefficients. Various 
formulae, often involving experimental input (ionization 
potentials and polarizabilities), have been proposed for 
this. 7 9 ' 12 However, this reliance on experimental data 
limited the set of elements that could be treated to those 
present in typical organic molecules. The third issue is 
that the Cq coefficients are kept constant during the cal- 
culation, and so effects of different chemical states of the 
atom or the influence of its environment are neglected. 
The fourth drawback, which we discuss later, is that the 
Cq/t 6 function diverges for small separations (small r) 
and this divergence must be removed. 

For a more widely applicable method a more consis- 
tent means of deriving the dispersion coefficients is re- 
quired. In 2006 Grimme published one such scheme, re- 
ferred to as DFT-D2. 17 In this approach the dispersion 
coefficients are calculated from a formula which couples 
ionization potentials and static polarizabilities of isolated 
atoms. Data for all elements up to Xe are available and 
this scheme is probably the most widely used method 
for accounting for dispersion in DFT at present. Whilst 
incredibly useful, DFT-D2 is also not free from prob- 
lems. In particular, for some elements arbitrary choices 
of dispersion coefficients still had to be made. For exam- 
ple, alkali and alkali earth atoms use coefficients that are 
averages of the previous noble gas and group III atom. 
Furthermore, the dispersion energy, E^ sp ^ is scaled ac- 
cording to the XC functional used and as a result the 
interaction energy of two well separated monomers is not 
constant but sensitive to the choice of XC functional. 

With the simple Cq/v q correction schemes the disper- 
sion correction diverges at short inter-atomic separations 
and so must be "damped". The damped dispersion cor- 
rection is typically given by a formula like 

£dis P = -£ /(t-ab, A, B)C 6 AB /ri B , (3) 

A,B 

where the damping function /(tab? A, B) is equal to one 
for large r and decreases i^disp to zero or to a constant 
for small r. We illustrate the divergence and a possible 
damping by the red curves in Fig. 5. How the damping is 
performed is a thorny issue because the shape of the un- 
derlying binding curve is sensitive to the XC functional 
used and so the damping functions must be adjusted so 
as to be compatible with each exchange-correlation or 
exchange functional. This fitting is also sensitive to the 
definition of atomic size (van der Waals radii are usu- 
ally used) and must be done carefully since the damp- 



ing function can actually affect the binding energies even 
more than the asymptotic Cq coefficients. 38 The fitting 
also effectively includes the effects of Cs/r 8 or Cio/r 10 
and higher contributions, although some methods treat 
them explicitly. 19,39 An interesting damping function has 
been proposed within the so-called DFT coupled cluster 
approach (DFT/CC), it has a general form that can ac- 
tually force the dispersion correction to be repulsive. 40 
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FIG. 5. Schematic illustration of a binding curve (i^tot) ob- 
tained from a dispersion corrected DFT calculation and the 
contributions to it from the regular DFT energy (.Edft) and 
the dispersion correction (i^disp)- The function — 1/r 6 used 
to model the dispersion correction diverges for small r (red 
dashed curve) and must be damped (solid red curve). Since 
the details of the damping strongly affect the position of the 
energy minima on the binding curve, the damping function 
needs to be fit to reference data. In addition, the damping 
function will be different for different XC functionals used to 
obtain Euft- 

Before leaving the simple Cq/t 6 schemes we note 
that although the approaches developed by Grimme are 
widely used, other parameterizations, which differ in, 
e.g., the combination rules used, are available. 40 43 Fur- 
thermore, functionals such as B97-D 17 and co>B97X-D 44 
have been specifically designed to be compatible with this 
level of dispersion correction. 



C. Step two — Environment-dependent Cq 
corrections 

A problem with the simple "DFT-D" schemes is that 
the dispersion coefficients are predetermined and con- 
stant quantities. Therefore the same coefficient will be 
assigned to an element no matter what its oxidation or 
hybridization state. However the errors introduced by 
this approximation can be large, e.g. the carbon Cq co- 
efficients can differ by almost 35% between the sp and 
sp 3 hybridized states. 9 Thus the emergence of methods 
where the Cq coefficients vary with the environment of 
the atom has been a very welcome development. We put 
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these methods on the second step of our stairway. They 
still use Equation 3 to obtain the dispersion correction 
and, as with the step 1 methods, some reference data 
(such as atomic polarizabilities) is used to obtain the Cq 
coefficients. 

We now discuss three step 2 methods: DFT-D3 
of Grimme, the approach of Tkatchenko and Schemer 
(vdW(TS)), and the Becke- Johnson model (BJ). The uni- 
fying concept of these methods is that the dispersion co- 
efficient of an atom in a molecule depends on the effective 
volume of the atom. When the atom is "squeezed", its 
electron cloud becomes less polarizable leading to a de- 
crease of the Cq coefficients. 

Grimme et al. 19 capture the environmental dependence 
of the Cq coefficients by considering the number of neigh- 
bors each atom has. When an atom has more neighbors 
it is thought of as getting squeezed and as a result the 
Cq coefficient decreases. This effect is accounted for by 
having a range of precalculated Cq coefficients for vari- 
ous pairs of elements in different reference (hybridization) 
states. In the calculation of the full system the appro- 
priate Cq coefficient is assigned to each pair of atoms ac- 
cording to the current number of neighbors. The function 
calculating the number of neighbors is defined in such a 
way that it continuously interpolates between the precal- 
culated reference values. Therefore if the hybridization 
state of an atom changes during a simulation, the Cq coef- 
ficient can also change continuously. Despite the simplic- 
ity of this approach, termed "DFT-D3", the dispersion 
coefficients it produces are pretty accurate. Specifically, 
based on the reference data of Meath and coworkers there 
is a mean absolute percentage deviation (MAPD) in the 
Cq coefficients of 8.4%. Moreover, the additional compu- 
tational cost is negligible since the number of neighbors 
can be obtained quickly. 

In 2009 Tkatchenko and Schemer 18 proposed a method 
which relies on reference atomic polarizabilities and ref- 
erence atomic Cq coefficients 45 to calculate the dispersion 
energy. These quantities are sufficient to obtain the Cq 
coefficient for a pair of unlike atoms. 46 To obtain environ- 
ment dependent dispersion coefficients effective atomic 
volumes are used. During the calculation on the system 
of interest the electron density of a molecule is divided 
between the individual atoms (the Hirshfeld partition- 
ing scheme is used) and for each atom its corresponding 
density is compared to the density of a free atom. This 
factor is then used to scale the Cq coefficient of a ref- 
erence atom which changes the value of the dispersion 
energy. The accuracy of the final (isotropic, averaged) 
Cq coefficients is quite high, with the MAPD on the data 
of Meath and coworkers being only 5.4%. However, it 
is not yet clear if the scaling of Cq coefficients with vol- 
ume will be accurate for more challenging cases such as 
different oxidation states in e.g. ionic materials. 

The most complex step 2 method is the BJ 
model, 13 16 ' 39 which exploits the fact that around an 
electron at r\ there will be a region of electron density 
depletion, the so-called XC hole. Even for symmetric 



atoms, this creates asymmetric electron density and thus 
non-zero dipole and higher-order electrostatic moments, 
which causes polarization in other atoms to an extent 
given by their polarizability. The result of this is an at- 
tractive interaction with the leading term being dipole- 
induced dipole. The biggest question in the BJ model is 
how to quantify the XC hole and in the method it is ap- 
proximated as the exchange-only hole (calculated using 
the Kohn-Sham orbitals). An average over the positions 
7*1 of the reference electron then gives the average square 
of the hole dipole moment, denoted (d%) which together 
with the polarizabilities a enters the formula for the dis- 
persion coefficient 

r AB = Q^b(4)a(4)b U] 

6 (d*) A OB + (d*) B a A ' [ } 

This is, in fact, very similar to the vdW(TS) formula 
but in vdW(TS) precalculated Cq coefficients are used 
instead of the hole dipole moments. The BJ model is very 
intriguing and several authors have studied how it relates 
to formulae derived from perturbation theory. 47 49 

In the BJ model the Cq coefficients are altered through 
two effects. First, the polarizabilities of atoms in 
molecules are scaled from their reference atom values ac- 
cording to their effective atomic volumes. Second, the 
dipole moments respond to the chemical environment 
through changes of the exchange hole, although this effect 
seems to be difficult to quantify precisely. The details of 
how to obtain atomic volumes and the exchange hole are 
known to affect the results obtained to some extent, 50 52 
but overall the accuracy of the asymptotic Cq coefficients 
obtained from the BJ model is quite satisfactory, with a 
MAPD of 12.2 % for the data of Meath and coworkers. 19 
Compared to the two previous step 2 approaches, how- 
ever, the computational cost is relatively high, in the 
same ballpark as the cost of a hybrid functional. 26 



D. Step three — Long-range density functionals 

The methods discussed so far require predetermined 
input parameters to calculate the dispersion interaction, 
either the Cq coefficients directly or the atomic polariz- 
abilities. Now we discuss approaches that do not rely on 
external input parameters but rather obtain the disper- 
sion interaction directly from the electron density. This, 
in principle, is a more general strategy and thus we place 
these methods on step 3 of our stairway. The methods 
have been termed non-local correlation functionals since 
they add non-local (i.e. long range) correlations to local 
or semi-local correlation functionals. 

The non-local correlation energy E^ 1 is calculated from 

Ef = JJ dr 1 dr 2 n(r 1 )(^(r 1 ,r 2 )n(r 2 ) . (5) 
This is a double space integral where n(r) is the electron 
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density and (p(ri, V2) is some integration kernel. The ker- 
nel is analogous to the classical Coulomb interaction ker- 
nel l/\r\ — r2 1 but with a more complicated formula used 
for ip(r\ 1 r2) with 0(l/\r\ — r2 1 6 ) asymptotic behavior. 
The formula has a pairwise form and thus ignores the 
medium between points r\ and r<i- Various forms for E£ l 
were proposed in the nineties 4 6 but were restricted to 
non-overlapping fragments. This rather severe limitation 
was removed by Dion et al in 2004 10 who proposed a 
functional form which can be evaluated for overlapping 
molecules and for arbitrary geometries. Within this ap- 
proach the XC energy E xc is calculated as 

E XC = E^ GA + E^ DA + Ef, (6) 

where the terms on the right hand side are the exchange 
energy in the revPBE approximation, 53 the LDA corre- 
lation energy, and the non-local correlation energy term, 
respectively. This method has been termed the van der 
Waals density functional (vdW-DF). vdW-DF is a very 
important conceptual development since it adds a de- 
scription of dispersion directly within a DFT functional 
and combines correlations of all ranges in a single for- 
mula. 

Since the development of the original vdW-DF a num- 
ber of follow up studies have aimed at understanding 
and improving the performance of the method. First, it 
turns out that vdW-DF tends to overestimate the long 
range dispersion interactions: Vydrov and van Voorhis 
have shown that when the Cq coefficients are evaluated 
the errors can be as large as ~30%. 54 To address this 
these authors proposed (computationally cheaper) func- 
tional that reduce the average errors by approximately 
50%. 54 57 The developers of the original vdW-DF tried 
to address its tendency to overbind at large separations 
by proposing vdW-DF2. 58,59 This functional, which in- 
volves changes to both the exchange and non-local corre- 
lation components tends to improve the description of the 
binding around energy minima, however, the Cq coeffi- 
cients it predicts are considerably underestimated. 56 Sec- 
ond, aside from the particular form of the choice of 
the exchange functional used in Equation 6 has received 
considerable attention. The original revPBE exchange 
functional chosen for vdW-DF often leads to too large 
intermolecular binding distances and inaccurate binding 
energies. To remedy these problems alternative "less re- 
pulsive" exchange functionals have been proposed 60-63 
which when incorporated within the vdW-DF scheme 
(Equation 6) lead to much improved accuracy. Of these 
the "optB88" and "optPBE" exchange functionals have 
been shown to offer very good performance for a wide 
range of systems. 60 ' 61 ' 64 ' 65 

Initially the non-local correlation functionals came, to 
a lesser or greater extent, 66 ' 67 with a higher computa- 
tional cost than GGAs or hybrid functionals. However, 
thanks to the work of Roman-Perez and Soler the com- 
putational cost of vdW-DF is now comparable to that of 
a GGA. 68 In addition, self-consistent versions of vdW- 



DF and several of its offspring are now implemented in 
widely distributed DFT codes such as Siesta, 68 VASP, 61 
QuantumESPRESSO, 69 and QChem. 66 Other ways to re- 
duce the cost of vdW-DF calculations have also been re- 
ported, for example, Silvestrelli has shown 70 ' 71 that uti- 
lizing Wannier functions to represent the electron den- 
sity allows for an analytic evaluation of the function- 
als proposed in the nineties. 5 ' 6 Particularly interesting 
is the local response dispersion (LRD) approach of Sato 
and Nakai 72 ' 73 which yields very accurate Cq coefficients. 
Overall, by not relying on external reference data, the 
step 3 approaches are less prone to fail for systems out- 
side of the reference or fitting space of step 2 methods. 
However, very precise calculations of the dispersion en- 
ergy are difficult and the formulae underlying vdW-DF 
and similar approaches can be somewhat complicated. 



E. Higher steps — Beyond pairwise additivity 

The main characteristic of the methods discussed so far 
is that they consider dispersion to be pairwise additive. 
As a consequence the interaction energy of two atoms or 
molecules remains constant no matter what material sep- 
arates them and all atoms interact "on their own" with 
no consideration made for collective excitations. 74 While 
such effects don't seem to be crucially important in the 
gas phase, especially for small molecules, they are impor- 
tant for adsorption or condensed matter systems where 
the bare interaction is screened. 75 We now briefly discuss 
some of the methods being developed which treat disper- 
sion beyond the pairwise approximation. The range of 
approaches is quite wide, from methods based on atom 
centered interactions, to methods involving density, to 
methods using electron orbitals. Because of the fresh- 
ness of most of the methods we discuss them together. 




FIG. 6. The long-range electron correlations of two isolated 
atoms can be described by an effective — 1/r 6 formula. In a 
condensed system the interaction is modified (screened) by 
the presence of the other electrons. 

Recently, the Axilrod-Teller-Muto 76 ' 77 formula has 
been used to extend the atom-centred pairwise ap- 
proaches to include three-body interactions. 19 ' 78 The 
triple-dipole interaction between three atoms A, B, and 
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C can be written as 



E 



ABC 



a 



ABC 3 cos(a) cos(/3) cos(7) + 1 



(7) 



' AB BC AC 



where a, /?, and 7 are the internal angles of the triangle 
formed by the atoms A, B, and C. The dispersion coef- 
ficient C$ BC is again obtained from reference data. No- 
tably, von Lilienfeld and Tkatchenko estimated the mag- 
nitude of the three body terms to be ~ —25% (destabilis- 
ing) of the two body term for two graphene layers. 78 Since 
the three-body interaction is repulsive for atoms forming 
an acute angled triangle, the interaction was found to 
be destabilising for stacked configurations of molecular 
clusters such as benzene dimers. 

The atom-centered approach can be used to approxi- 
mate the many-body dispersion interaction using a model 
of coupled dipoles. 79 ' 80 Here quantum harmonic oscilla- 
tors with characteristic frequencies occupy each atomic 
position and the dispersion interaction is obtained from 
the shifts of the frequencies of the oscillators upon 
switching on their interaction. The model has been ap- 
plied outside DFT for some time 81 85 and has recently 
been used to show the non-additivity of dispersion in 
anisotropic materials. 86 The initial applications of this 
approach termed many-body dispersion (MBD) are quite 
promising and suggest that the higher-order dispersion 
terms can be important even for systems such as solid 
benzene. 80 However, getting accurate relations between 
atoms and oscillator models seems to be rather challeng- 
ing, especially when describing both localized and delo- 
calized electrons. 

In the context of DFT, orbitals can be used to calcu- 
late the correlation energy using the adiabatic-connection 
fluctuation-dissipation theorem (ACFDT). 6 ' 87 The par- 
ticular approach which has received the most attention 
recently is the so-called random phase approximation 
(RPA). 88 ' 89 Results from RPA have been very encour- 
aging and it exhibits, for example, a consistent accu- 
racy for solids 90 ' 91 and the correct asymptotic descrip- 
tion for the expansion of graphite, a feature which the 
pairwise methods fail to capture. 37 ' 92 RPA is analogous 
to post-HF methods and indeed direct links have been 
established. 93 ' 94 Unfortunately, it also shares with the 
post-HF methods a high computational cost (the cost in- 
creases approximately with the fourth power of the sys- 
tem size) 95 ' 96 and slow convergence with respect to basis 
set size. 97 101 Approaches going beyond RPA are also re- 
ceiving attention, for example the second-order screened 
exchange 102 or single excitations corrections. 103 

Another method involving orbitals to obtain the corre- 
lation energy is to combine DFT with post-HF methods 
in the hope of exploiting the benefits of each approach. 
While the post-HF methods can describe long range in- 
teractions accurately (albeit expensively) using orbitals, 
DFT is efficient for the effective description of the short 
range part of the interactions. So called range separated 
methods are an example of this approach where long 



range correlations are treated by, e.g., second-order per- 
turbation theory, the coupled cluster method, or RPA. 
Range separated functionals can be very accurate, can 
account for many-body interactions, and can deliver rel- 
atively fast convergence of the correlation energy with re- 
spect to basis set size. 104 109 Another example from this 
class are the so called double hybrid functionals, which 
include Fock exchange and a second order perturbation 
theory correlation contribution. Generally the coeffi- 
cients for these contributions are fitted to reproduce ref- 
erence data. 110 ' 111 Since the original double hybrid func- 
tionals were more concerned with reaction energies and 
barrier heights than dispersion, dispersion interactions 
were actually underestimated. Newer functionals, e.g., 
mPW2PLYP-D, add a dispersion correction in the sense 
of step 1 or 2 methods, however, they are computation- 
ally very expensive and at present prohibitively expensive 
for most condensed phase and surface studies. 33 ' 111 



F. Summary 



The higher steps of the stairway with the promising 
schemes that go beyond the pairwise approximation close 
our overview and classification of methods for treating 
dispersion within DFT. In Table I we summarize some 
of the key aspects of the most relevant schemes, such as 
what properties they use to obtain the dispersion con- 
tribution (via, e.g., the Cq) and the approximate relative 
computational cost. The computational cost of a method 
is, of course, a key factor since no matter how accurate 
it is, it will not be used if the computational cost is pro- 
hibitive. In this regard the cost of the simple pairwise 
corrections on step 1 is essentially zero compared to the 
cost of the underlying DFT calculation and these meth- 
ods can be recommended as a good starting point for 
accounting for dispersion. Because of the correct asymp- 
totic behavior (at least for gas phase molecules) 37 they 
are preferable to methods from ground. Compared to the 
step 1 methods the DFT-D3 and vdW(TS) schemes don't 
add a significant computational cost, whereas the vdW- 
DF approach increases the computational time by about 
50% compared to a GGA calculation. Since the accuracy 
of these approaches tends to be higher than that of the 
step 1 methods, they should be preferred over the step 1 
corrections. However, when applied to a particular prob- 
lem, the accuracy of any method used should be tested 
or verified against experimental or other reference data 
since any approach can, in principle, fail for a specific 
system. The methods on higher steps are mostly in de- 
velopment and currently the ACFDT or range separated 
functionals require a computational cost two to four or- 
ders of magnitude higher than a GGA calculation which 
limits their applicability. 



8 



TABLE I. Summary of some of the most widely used methods for capturing dispersion interactions in DFT, ordered according 
to the "Step" they sit on in the stairway to heaven (Fig. 4). Information on the reference used for the C6 coefficients, what the 
C6 depend on, and the additional computational cost of each approach compared to a regular GGA calculation is reported. 



Method 


Step 


Reference for C6 


C6 depend on 


Additional 
computational cost a 


Ref. 


Minnesota 





None 


N/A 


None 


e.g. 34 


DCACP 





None 


N/A 


Small 


11 


DFT-D 


1 


Various 


Constant 


Small 


e.g. 110 


DFT-D3 


2 


TDDFT 


Structure 


Small 


19 


vdW(TS) 


2 


Polarizabilities 
and atomic Cq 


Atomic volume 


Small 


18 


BJ 


2 


Polarizabilities 


Atomic volume, X hole 


Large 


14 


LRD 


3 


Cq calculated 


Density 


Small 


72 


vdW-DF 


3 


Cq calculated 


Density 


^50% 


10 


Dbl. hybrids 


4 


None or as "-D" 


Orbitals 


Large 


110 



a The BJ model and the double hybrids (labelled "Dbl. hybrids") 
are more computationally expensive than the simpler "DFT-D" 
methods 26 and the calculation of the correlation energy in 
vdW-DF leads to a ^50% slow-down compared to a GGA 
calculation when done efficiently. 68 
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III. GENERAL PERFORMANCE 

We have commented in passing on how some of the 
methods perform, mainly with regard to the dispersion 
coefficients. Now we discuss accuracy in more detail by 
focussing on binding energies for a few representative gas 
phase clusters, a molecular crystal, and an adsorption 
problem. Before doing so, it is important to emphasise 
two points. First, establishing the accuracy of a method 
is not as straightforward as it might seem since it must 
be tested against accurate reference data (e.g. bind- 
ing energies). Experimental data can be inaccurate and 
is rarely directly comparable to theory since quantum 
nuclear and/or thermal effects, which affect the experi- 
mental values, are often neglected in simulation studies. 
Theoretical reference data, on the other hand, is more 
appropriate but rather hard to come by since accurate 
reference methods such as post-HF methods (MP2, CC) 
or quantum Monte Carlo (QMC) are so computationally 
expensive that they can generally be applied to only very 
small systems (tens of atoms). 30 Second, the results ob- 
tained with many dispersion based DFT approaches are 
often strongly affected by the fitting procedure used when 
combining the long range dispersion interaction with the 
underlying exchange or XC functionals. 



A. Gas phase clusters 

Because of growing computer power and algorithmic 
improvements, systems of biological importance such as 
DNA and peptides can now be treated with DFT. As 
discussed in the introduction, dispersion plays an impor- 
tant role in stabilising the folded state and thus disper- 
sion corrected functionals are required. However, since 
accurate reference data can only be obtained for systems 
with tens of atoms, it's just not possible to obtain refer- 
ence data for DNA or a protein directly and fragments of 
the large molecules are used to build test sets of binding 
energies and geometries instead. One such test set 112 114 
is the popular S22 data set of Jurecka et al.? 2 It is use- 
ful since it contains 22 different dimers with a range of 
weak bonding types and with a wide range of interac- 
tion energies. Results for some of the methods are re- 
ported in Table II; specifically we report mean absolute 
deviations (MAD) and mean absolute percentage devi- 
ations (MAPD). Included in Table II are results from 
LDA and semi-local PBE. These functionals fail to re- 
produce the correct interaction energies, yielding MADs 
of ~10 kJ/mol and MAPDs in excess of 30%. Much im- 
proved performance is observed with even the simplest 
pairwise corrections schemes. Indeed, on all steps of the 
stairway above ground there is at least one method with 
a MAD below 1.5 kJ/mol (or 0.4 kcal/mol, ~16 meV) 
and MAPD close to 5%. This is very good agreement 
with the reference data and a clear indication of the re- 
cent improvements made in the DFT-based description of 
dispersion. As can be seen, the results obtained on lower 



steps of the stairway can surpass those from methods on 
higher steps. This is not that surprising since some of the 
methods reported in Table II were actually developed by 
fitting to the S22 data set itself. 



TABLE II. Mean absolute deviations and mean absolute per- 
centage deviations of different dispersion based DFT meth- 
ods on the S22 data set of Jurecka et al. 22,115 Results are 
in kJ/mol for MAD and percent for MAPD. Results from 
second order M0ller-Plesset perturbation theory (MP2) - a 
widely used post-HF method - are shown for comparison as 
are results from the LDA and PBE functionals as examples 
of methods that don't treat dispersion explicitly. 



Method 


MAD 


MAPD 


Ref. 


Ground 








LDA 


8.66 


30.5 


116 


PBE 


11.19 


57.6 


116 


M06-2X 


2.52 


11.3 


117 


LAP 


2.51 


7.8 


35 


Step 1 








B97-D2 


1.51 


7.3 


44 


o;B97X-D 


0.76 


5.6 


44 


Step 2 








BLYP-D3 


0.96* 




19 


PBE-vdW(TS) 


1.25** 


9.2** 


18 


rPW86-BJ 


1.50 


6.1 


51 


Step 3 








LC-BOP+LRD 


0.86 


4.6 


73 


vdW-DF 


6.10 


22.0 


67 


optB88-vdW 


1.18 


5.7 


60 


vdW-DF2 


3.94 


14.7 


118 


VV10 


1.35 


4.5 


118 


Higher steps 








B2PLYP-D3 


1.21* 




19 


cjB97X-2 


1.17 


8.8 


119 


PBE-MB 




5.4** 


80 


RPA 


3.29** 




101 



post-HF 

MP2 3.58 19.4 115 

* using the original reference values of Jurecka et al. 
(Ref. 22) 

** using the values of Takatani et al. (Ref. 120) which 
are very similar to those of Podeszwa et al. (Ref. 115) 

Small water clusters, in particular water hexamers, 
are an interesting system where dispersion plays a sig- 
nificant role (see, e.g. Refs. 121-125). The four 
relevant isomers (known as "book", "cage", "cyclic", 
"prism", Fig. 7) all have total energies that differ by 
<1.3 kJ/mol per molecule according to accurate post- 
HF methods. 60 ' 121 While the "prism" and "cage" iso- 
mers are preferred by post-HF methods, standard XC 
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Cage 
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1 0-5 



1-0.5 



Book 



Cyclic 



« -l 



X Reference 
+ revPBF-vdW 
• optB88-vdW 
■ BLYP-vdW(Silv.) 
A PBE-vdW(TS) 
% PBE 



Prism Cage Book Cyclic 



FIG. 7. The four low energy isomers of the water hexamer 
are an important example where dispersion interactions play 
a significant role. In the panel on the right binding energy dif- 
ferences per water molecule as predicted by different XC func- 
tional are compared to accurate reference data. 60 ' 121 While 
the reference method predicts the "prism" and "cage" isomers 
to be more stable, standard XC functionals such as PBE give 
the "book" and "cyclic" clusters as the ones with lower energy. 
It was shown in Ref. 121 that the agreement with the refer- 
ence calculations is improved when dispersion interactions are 
accounted for, shown here by results with various dispersion 
correction schemes. 



functionals find the more open "cyclic" and "book" iso- 
mers to be more stable. Recent work reveals that dis- 
persion impacts profoundly on the relative energies of 
the isomers and improved relative energies are obtained 
when dispersion is accounted for. Indeed the water 
hexamers have become an important test system with 
techniques such as DFT-D, 121 vdW(TS), 121 vdW-DF, 126 
optB88-vdW, 60 and the modified vdW-DF approach of 
Silvestrelli (BLYP-vdW(Silv.)) 127 all having been ap- 
plied. The performance of these schemes in predicting 
the relative energies of the water hexamers is shown in 
Fig. 7, where one can also see the improvement over 
a standard functional such as PBE. One general point 
to note, however, is that even without a dispersion cor- 
rection certain functionals can already give quite accu- 
rate absolute binding energies for systems held together 
mainly by hydrogen bonding and adding dispersion cor- 
rections can actually lead to too large absolute binding 
energies. This is indeed the case for the water hexam- 
ers where, for example, the vdW(TS) correction to PBE 
gives binding energies for the hexamers that are too large 
by ~4 kJ/mol per molecule. 



B. Molecular crystals — solid benzene 

Molecular crystal polymorph prediction is another im- 
portant area where dispersion forces can play a criti- 
cal role, and even for molecular crystals comprised of 
small molecules several polymorphs often exist within a 
small energy window. 128 ' 129 Identifying the correct ener- 
getic ordering of the polymorphs can therefore be a strin- 
gent test for any method. A large number of molecular 
crystals have been characterised experimentally (see, e.g. 



Ref. 130) and this, therefore, could serve as a rich testing 
ground for DFT-based dispersion schemes. 

Solid benzene is one of the most widely examined test 
systems, with studies focussing on the experimental den- 
sity and the cohesive energy of the crystal. 131 The lat- 
ter is obtained from the experimental sublimation en- 
ergy from which the effects of temperature and quan- 
tum nuclear effects must be subtracted giving a value 
in the 50 to 54 kJ/mol per molecule range. 132 Calcula- 
tions with PBE give an abysmal cohesive energy of only 
10 kJ/mol per molecule and a volume ^30% too large 
compared to experiment. 133 This huge difference between 
theory and experiment suggests that dispersion is impor- 
tant to the binding of the crystal and indeed the error 
is greatly reduced when even rather simple schemes are 
used. For example, PBE-D2 and the DCACP poten- 
tials give estimates of the lattice energy (55.7 kJ/mol 133 
and 50.6 kJ/mol, 36 respectively) and densities that are 
within 10% of experiment. While the PBE-vdW(TS) 
scheme overbinds slightly (66.6 kJ/mol), it has been sug- 
gested that the many-body interactions are quite im- 
portant for this system and their inclusion reduces the 
PBE-vdW(TS) value by 12 kJ/mol. 80 Non-local vdW- 
DF overestimates the cell volume by ~10% but gives 
a rather good cohesive energy of 58.3 kJ/mol. 116 vdW- 
DF2 reduces both the binding energy (to 55.3 kJ/mol) 
and the error in the cell volume, which turns out to be 
~3% larger than the experimental value. 116 Using func- 
tionals proposed to improve upon the overly repulsive 
behavior of vdW-DF, such as optB88-vdW leads to an 
improvement of the reference volume (~3% smaller than 
experiment), but the cohesive energy is overestimated, 
being 69.4 kJ/mol. 116 Again, the overestimated cohe- 
sive energy may result from missing many-body inter- 
actions. The RPA has also been applied to this system 
and while the density is in very good agreement with ex- 
periment, the cohesive energy is slightly underestimated 
at 47 kJ/mol. 132 Overall, the improvement over a semi- 
local functional such as PBE is clear but more applica- 
tions and tests on a wider range of systems are needed 
to help establish the accuracy of the methods. 



C. Adsorption — benzene on Cu(lll) 

Adsorption on solid surfaces is another area where 
great strides forward have recently been made with re- 
gard to the role of dispersion. Dispersion is of obvious 
importance to the binding of noble gases to surfaces but it 
can also be important to chemisorption 134 and, e.g., wa- 
ter adsorption. 64 Indeed for weakly chemisorbed systems 
dispersion forces attain an increased relative importance 
and one can estimate from various studies (e.g., Refs. 135 
and 136) that dispersion contributes about 4 to 7 kJ/mol 
to the adsorption energy of a carbon-sized atom; not a 
negligible contribution. 

Of the many interesting classes of adsorption system, 
organic molecules on metals have become a hot area of re- 



11 




FIG. 8. In molecular crystals the molecules can be packed in 
a range of orientations and at a range of distances from each 
other. This makes molecular crystals a more challenging class 
of system than gas phase clusters, and, indeed, even for the 
relatively simple case of benzene shown here, few methods are 
able to very accurately describe the cohesive properties of the 
crystal. The benzene molecules in the central unit cell are 
shown with large balls, all other benzenes in the neighboring 
unit cells are shown in wireframe. 



search and for such systems dispersion must be included 
if reasonable adsorption energies and structures are to 
be obtained. Benzene on Cu(lll) is an archetypal widely 
studied system for organic adsorbates on metals. It is also 
an interesting system because it illustrates how difficult 
it is sometimes to use experimental reference data, which 
for this system has proved to be somewhat of a "moving 
target". Indirect estimates obtained from temperature 
programmed desorption initially placed the adsorption 
energy at 57 kJ/mol, 137 however, a recent reinterpreta- 
tion of the experiment moved the adsorption energy up to 
the 66 to 78 kJ/mol range. 138 When looking at this sys- 
tem with DFT, not unexpectedly, PBE gives very little 
binding (5 kJ/mol). 139 Most of the step 1 and 2 meth- 
ods overestimate the binding, 140 for example, PBE-D2 
and PBE-vdW(TS) give 97 and 101 kJ/mol, respectively. 
The exception is the work of Tonigold and GroB 135 where 
a value of 59 kJ/mol was obtained, based on dispersion 
coefficients obtained by fitting to post-HF data of small 
clusters. Recently, Ruiz et al. 138 approximately included 
the many-body effects in calculating the Cq coefficients 
in the vdW(TS) scheme, which reduced the predicted ad- 
sorption energy to 88 kJ/mol. vdW-DF gives adsorption 
energies of about ~53 kJ/mol which underestimates both 
the old and new reference data. 139 ' 141 It is clear from this 
and other systems that adsorption is very challenging for 
dispersion-based DFT methods at present. The "DFT- 
D" methods face the problem of obtaining the Cq coef- 
ficients for the atoms within the surface of the solid and 
the pairwise methods neglect many-body effects. Both of 
these issues will require much more consideration in the 
future. 



IV. FINAL REMARKS 

An enormous amount of progress has been made over 
the last few years with the treatment of dispersion forces 
within DFT. It is now one of the most exciting and thriv- 
ing areas of development in modern computational ma- 
terials science and an array of methods has been devel- 
oped. Here we have introduced and classified some of the 
main, often complementary, approaches. We have also 
discussed how some schemes perform on a selection of 
systems, including gas phase clusters, a molecular solid 
and an adsorption problem. These examples, and the 
many others in the literature, demonstrate that the range 
of systems which can now be treated with confidence with 
DFT has been greatly extended. Connected with this, 
the variety of systems and materials to which dispersion 
forces are now thought to be relevant has grown substan- 
tially. As a result the mantra that "dispersion forces are 
not important" is heard less often now and there is much 
less of a tendency to sweep dispersion forces under the 
rug. 

We have seen how many methods have been applied 
and tested on gas phase molecular clusters, for which 
there are now several approaches that can yield very 
high accuracy. However, a key contemporary challenge 
is the need to develop methods that will be accurate for 
both gas phase molecular systems and problems involv- 
ing condensed matter such as adsorption. Here, much 
work remains to be done with regard to the develop- 
ment of DFT-based dispersion techniques for condensed 
matter as well as in simply better understanding how 
current techniques perform in e.g. adsorption systems 
where there may be strong polarization effects. In this 
regard, approaches based on the ACFDT such as RPA 
look promising. However, given their high computational 
cost and complex set-up it is likely that for the foresee- 
able future such methods will mainly be useful for tour de 
force reference style calculations rather than for routine 
studies. 142 

Looking to the future, the efficient description of many- 
body correlation effects in metals and other solids is an 
important unresolved issue. Even simply better under- 
standing their importance for different systems would be 
useful. For example, although many body correlation ef- 
fects should be important in solids, the vdW-DF method 
which neglects them can perform surprisingly well for 
solids. 61 Another problem closely related to dispersion is 
the issue of screening which differs substantially in solids 
and molecules. For molecules, approaches that neglect 
screening of exchange, e.g., the so-called long-range cor- 
rected (LC) exchange functionals 143 are beneficial. In- 
deed, the LC functionals improve many properties of 
molecules such as electrostatic moments which in turn 
decreases an important source of errors in vdW bonded 
systems. 55 ' 119 ' 144 In solids, especially metals or semicon- 
ductors, the interaction of electrons is significantly mod- 
ified by the presence of the other electrons and it is 
essential to capture this XC effect. Indeed, no screen- 
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ing of the Fock exchange leads, for example, to over- 
estimated band gaps of semiconductors. 145 The effect is 
also significant for correlation, for example, when two ho- 
mogeneous electron gas spheres are taken from vacuum 
into a homogeneous electron gas background their dis- 
persion interaction is reduced to about a fifth of its orig- 
inal value. 146 ' 147 Since screening is system dependant it 
is unclear how to treat solids and molecules on the same 
footing in a computationally efficient manner. A fur- 
ther issue, which is often of minor importance but should 
be accounted for when high accuracy is being sought is 
anisotropy of the dispersion coefficients. Isotropic dis- 
persion coefficients seem to be a good first approxima- 
tion since the anisotropy for molecules is on the order 
of 10%. 148 Anisotropy can, however, become an issue 
for highly anisotropic and polarizable objects 86 and to 
establish the overall importance would be useful. The 
anisotropy can be included in some schemes, as has been 
done, for example, in the BJ model and the vdW(TS) 
approach. 80 ' 149 

Many factors have prompted the recent progress with 
DFT, with one key factor being the parallel development 
of post-HF methods. This has provided the accurate ref- 
erence data which has served to both shine light on prob- 
lems with existing XC functionals and against which new 
methods can be proved. The fact that some of this refer- 
ence data has been easily accessible - such as the S22 data 
set 150 - has also helped. However, as stressed above, an 
important challenge nowadays is to develop methods that 



are accurate for solids and for adsorption. Unfortunately 
for these systems accurate reference data are scarce and, 
indeed, urgently needed either from experiment (e.g. 
microcalorimetry 151 for adsorption) or higher level elec- 
tronic structure theories (e.g. post-HF methods, 152 157 
approximations of ACFDT, QMC). It is encouraging that 
progress is being made in both of these areas. It is 
also very encouraging that condensed phase reference sys- 
tems are beginning to emerge, such as ice, LiH, water on 
LiH and water on graphene. 63 ' 154 ' 155 ' 157 164 By tackling 
these and other reference systems with the widest pos- 
sible range of techniques we will better understand the 
limitations of existing dispersion-based DFT approaches, 
which will aid the development of more efficient and more 
accurate methods for the simulation of materials in gen- 
eral. 
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